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Abstract 

It is well known that the Poisson-Boltzmann (PB) equation yields the exact coun- 
terion density around charged objects in the weak coupling limit. In this paper we 
generalize the PB approach to account for coupling of arbitrary strength by making 
use of Tsallis q-exponential distributions. Both the weak coupling and the strong 
coupling limits are reproduced. For arbitrary coupling we also provide simple ana- 
lytical expressions which are compared to recent Monte Carlo simulations by A. G. 
Moreira and R. R. Netz [Europhys. Lett. 52 (2000) 705]. Excellent agreement with 
these is obtained. 
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1 Introduction 



Many biopolymers, as DNA and actin, are strongly charged and the expla- 
nation of the screening of their charges by counterions in solution is of large 



relevance (for a recent review, see 



l||). A detailed knowledge of these screening 



effects not only allows to determine properties of individual polyelectrolytes, 
but also forces acting between them and the thermodynamic and transport 
properties of their solutions 0,0,00, 5|. The traditional Poisson-Boltzmann 
(PB) equation obtained from Boltzmann-Gibbs (BG) thermostatistics has long 
been known to fail to explain the concentration profiles of multivalent counte- 
rions close to a charged macroion 2j,|6j]. The Perel-Shklovskii CPS) theory 2\ 
based in the strongly correlated liquid and the Wigner crystal pj concepts has 
provided a qualitative explanation of this failure and also some quantitative 
results that have been recently confirmed by numerical simulations a,l2|. The 
approach has been succesfully applied to explain charge inversion. However, 
an analytical expression for the counterion concentration profile valid for the 
whole spatial range and for any ionic number Z has not yet been proposed. The 



3 



recent theory of Moreira and Netz [lOj] making use of field-theory methods, 
rigorously establishes the strong coupling limit for the counterion distribution 
besides a charged planar wall. The expression provided in this limit by the 
authors agrees with their Monte Carlo simulations, which are also carried for 
moderate coupling [10]. For the weak coupling regime there exists also den- 



12 



sity functional approaches [11 1 and methods involving integral equations 
that provide accurate results. However, a systematic and unified theory is still 
lacking for coupling of arbitrary strength. 



In the present paper we attempt to provide such a theory. The theory exactly 
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reproduces both the weak and strong coupling limits being also consistent with 
classical rigorous results. One of the main advantages of the theory is that 
strikingly simple expressions are provided for arbitrary coupling in excellent 



agreement with Monte Carlo simulations by Moreira and Netz [lOj . The PB 
equation is generalized to account for the correlations making use of the Tsallis 
q-exponential distributions. 

We present first results of the classical PB theory. Let us consider a massive 
insulating macro ion having a negative surface charge density — a. Assume 
further that the only counterions in solution are those dissociated from the 
surface so that the system macroion+solution is neutral. When there are no 
charge correlations (monovalent counterions and low a) the mean field PB 
equation 

Aip = -AnlN(0)e-^ (1) 

describes accurately the counterion concentration profile that screens the charge 
of the macroion. In Eq.(l) ip = —/3eZ[<f)(x) — 0(0)] {(3 = 1/kT is the room La- 
grange parameter, e is the electron charge and Z is the ionic number) N(0) and 
0(0) being the counterion concentration and the electric potential at the sur- 
face of a charged object respectively. In Eq. (1), I = Z 2 Ib, where Ib = e 2 / DkT 
is the Bjerrum length, which is a characteristic threshold for repulsion forces. 
The relative permittivity for water is D w 80. Eq.(l) results from introducing 
in the Poisson equation Aip = —4nlN(x) the Boltzmann equation 

N{x) = N{0)e~^ (2) 

which relates the counterion concentration at a position x to that in the vicin- 
ity of the surface, considering a planar geometry. Two boundary conditions 
are to be supplemented to Eq.(l) 
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dip _ 2 
dx x=0 A 

dx 

x— >oo 

where A = D / {2nf3aZe) is the Gouy-Chapman length 



(3) 
(4) 

1]. Eqs.(l), (3) and (4) 



completely specify the boundary value problem. Its solution is 



N(x) 



1 1 



(5) 



2ttI (x + A) 2 

This profile correctly describes the counterion density in the weak coupling 
limit. Out of this limit, significant departures (of importance for many phe- 
nomena) occur. It has been shown recently that in the limit of strong coupling, 
the counterion concentration profile decays exponentially with the distance to 
the macroion surface 



N(x) 



(6) 



The mean field PB approach cannot account for this profile. Furthermore, 
although extensive Monte Carlo calculations have been performed for this 



system 



10l | . a theory for arbitrary coupling is still lacking. 



2 Generalized mean field theory 

The exact evaluation of the Gibbs canonical partition function for the many 
(interacting) particle system provides the best description of the counterion 
distributions. However, an analytical closed form for this canonical partition 
function is inviable and requires computational resources in most cases. The 
PB theory above outlined is a valid approximation when correlations are not 
important and provides quite useful expressions strictly valid only in the weak 
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coupling limit. 



We can wonder if a mean field description can still be possible when coupling 
of arbitrary strength is incorporated. If one assumes a mean field potential 
to hold, its value must then be allowed to fluctuate. We can assume that 
fluctuations of the mean field value of the electrostatic potential are distributed 
according to a \ 2 ( or gamma) probability distribution function of order n, 

fW) 

n 

Then, all possible values for the mean field value of the electrostatic potential 
must be weighted by this distribution and the canonical partition function 



becomes 



13] 



f^')e-^'d^={l-(l- q )^/py- 



which coincides with the Tsallis q-exponential distribution (here q = 1 + 2/n). 
These expressions hold rigorously for q > 1 although for q < 1 there exists 



an analogous treatment 13j]. It is to be noted, however, than in previous 
considerations the fluctuating variable is the Lagrange parameter and not the 
energy. We consider here fluctuations in energy in a rather formal way due 
to the fact that energy is still approximated to a mean field value beyond its 
traditional domain. Fluctuations of arbitrary strenght must hence be included. 
j3 is proportional to the room Lagrange parameter (3 (i.e. the Boltzmann-Gibbs 
Lagrange parameter) and related through a constant k q which can depend on 
q. Clearly, if q = 1 in this expression then f(ip') = and the Boltzmann 

factor corresponding to the PB mean field theory is regained, q is the Tsallis 
entropic parameter, which is related in the present problem to the coupling 
constant entering in the Hamiltonian. 
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For coupling of arbitrary strength Eq.(2) can then be generalized in the fol- 
lowing way 

N(x) = N(0)e-^ /l3 = N(0) [l - (1 - q)]fy/p\ ^ (9) 

In this work q is restricted to the range < q < 1. When q = 1, the Boltzmann 
equation, Eq. (2) is regained. 



q-exponentia 
applications 



distributions are extremely interesting and have proven many 



14j . These can be obtained directly by maximizing the Tsallis 



entropic form S q 



15] 



T a v q - 1 

S q = k^^- (10) 



16] 



imposing the following constraint for the biased average of the energy 

Y n v q e ■ 

The q-exponential distributions have been applied succesfully to systems in 
which some degree of fractality is present, and are used here tentatively pro- 
vided that correlated structures possess fractal features as it is well known, 
for example, from studies of critical phase transitions. 

It is important to note that it is not the validity of BG thermostatistics what 
is at issue here. As indicated above, BG thermostatistics provides the best 
description of the physics involved. What we suggest is to generalize the PB 
mean field theory, which is based on the factorization of the many particle 
probability distributions. When correlations are important, this factorization 
cannot be carried out in the framework of BG thermostatistics. It is proposed 
that the factorization in terms of single particle probability distributions is, 
indeed, still possible if one uses the Tsallis generalized framework, considering 
the appropriate probability distributions. Since a mean value for the electro- 
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static potential is still assumed to hold, it must be allowed to strongly fluctu- 
ate and this is what we have done above by assuming a chi-square probability 
distribution for the fluctuations. 



Tsallis statistics has proven a valuable tool in dealing with long range interac- 
tions. The existence of screening in the case of Coulomb systems and, hence, 
of a thermodynamic limit Q| is not, in principle, a matter for worrying about 
concerning its application. It is known, for example, that the Tsallis canon- 
ical partition function can be derived assuming that a certain hierarchy of 
structures (for example a polydisperse distribution of particle sizes) is present 



with a given probability |18(. The interaction between these structures, con- 
sidered at the ensemble level, has a thermodynamic limit. In our system, the 
structures could be understood as those composed by each counterion and a 
number of shells of its screening atmosphere. The mean field theory here is 
not applied directly to each counterion moving in a potential created by the 
others, but to a structure containing the counterion and a portion of screen- 
ing atmosphere of averaged size (which is mediated by the surface charge and 
the ionic valence and is contained in the q-dependence). In this framework 
q pa 1 accounts for disordered screening shells corresponding to low charged 
ions (and low surface charge concentration) whereas q <C 1 describes highly 
charged charged systems with screening shells having some degree of order. 
These screening shells are of fractal character. Order enters in the present 
description by remembering that the parameter q tends to privilege certain 
microstates again others between all available ones. Interestingly, a strongly 
correlated liquid has been suggested to form a Wigner crystal at the macroion 
surface [3|, implying favoring certain concrete microstates against all those 
allowed in the disordered bulk phase. 
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In the theory developed here the Tsallis canonical partition function is hence to 
be considered as the correct way of factorizing the many-particle distribution 
of the exact BG canonical partition function in single-particle (mean field) 
contributions when the pair interaction term is considered. The q-dependence 
and the structure of the Tsallis canonical partition function corresponds then 
to a more general averaging which includes correlations in a natural way. 
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It is also worth noting that there exists also an extensive approach 19] that 
leads to Tsallis q-exponential distributions. This is of interest to establish a 



variational principle following analogous lines to those of Ref. 
the following generalized PB equation 



3- 



It leads to 



A?/; = -4vr/iV(0)e- w//3 = -AnlN{0) [l - (1 - q)j3if>/p] 1-9 (12) 
which can also be easily obtained by replacing Eq. (9) in the Poisson equation. 

The boundary value problem specified by Eqs.(3), (4) and (12) can now be 
analytically solved in an analogous way to that followed in the PB approach, 
which led to Eq. (5). The counterion concentration profile is then given by 

M\ (2-g)^ , qPx\~ Vq n ^ 



As can be observed, for x = 0, the counterion concentration is 



"<°> = mt4 (14) 



An exact result obtained with the classical theory and which must also hold 



here comes from the contact theorem 



2l| which states that JV(0) = l/(2vr/A 2 ). 



This allows to know (3 in terms of the actual Lagrange parameter (3 thus 
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removing the arbitrariness of this parameter. We obtain 



P 



2-q 



(15) 



and, therefore, the counterion concentration profile can be finally written as 

-2/q 



N(x) 



1 f i+ qx 



(16) 



2tt/A 2 V (2-q)\J 
As can be seen, for q — 1 (weak-coupling limit) Eq.(5) is regained. The coun- 



terion concentration profile is normalized independently of q 



oo la 
N(x)dx = — — = — - 



(17) 



In the strong coupling case (q — > 0) Eq. (16) describes a pure exponential 
decay 



N{x) 



1 



2^P eXP V A 



x 



;is) 



in agreement with the strong-coupling theory by Moreira and Netz, Eq. (6). 

Since the whole system is neutral, the apparent surface charge density at 
position x can be calculated by making use of Eq.(16) as 



a* = Ze 



roc 

/ N{x')dx' 

Jx 



a 



1 + 



qx 



(2-q)\ 



(19) 



It is worthy noting that for strong coupling with g/2 < 1 and close to the 
macroion surface x/X 2/q the counterion concentration profile, Eq. (16) 
can be approximated by 



N{x) 



2x 



exp 



exp 



A V 2 



(20) 



2vr/A 2 V (2 - q)X J 2vr/A 2 

We see that the strong coupling solution is the leading term of this expan- 
sion. This behavior at low distance from the surface is totally consistent with 
previous heuristic considerations I22 1 . At large distances from the macroion 
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surface, x = xq + x (being i <C i , a;o large) we have approximately 

N(x) « N(x ) (J^j « (l + ~) ' (21) 

with A = Xoq. PB behavior is thus regained if the surface is taken as if it 
were at Xq. When lowering q, the position xq must increase in order to be 
valid the expansion outlined. In the limit q — > 0, xn goes to infinity. This 



asymptotic behavior is consistent with PS theory 



221 ] in which the boundary 



condition is displaced. There is no previous analytical theory, however, capable 
of reproducing the crossover between strong coupling and PB like behaviors 
1(1 22^ . We achieve this here, through Eq. (16), for arbitrary coupling strength 
and position. 



3 Comparison with Monte Carlo simulations 



Moreira and Netz |10( have performed extensive Monte Carlo calculations on 
the system under consideration. These are based on the specification of the 
coupling constant S, which enters in the Hamiltonian of their BG canonical 
partition function. The value of this constant in terms of the system parame- 
ters is l^J 

5 = 2nZ%a (22) 

This coupling constant gives the strength of the coupling and depends on 
the surface charge density of the macroion and on the ionic number Z of the 
counterions. 

In Fig.l we plot the counterion concentration profiles for the values of q indi- 
cated in the figure. Good fittings are obtained to Monte Carlo simulations by 
the authors of Ref. for different coupling strengths and also for moderate 
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coupling (q = 0.2) where there exists no previous systematic theory to the best 
of our knowledge. The entropic parameter is directly related to the coupling 
constant considered in the simulations. In Table 1 we indicate the values for 
the parameter q related to the values of H in the simulations. 

The data in Table 1 can be fitted accurately to a function q = 1/(1 + aH) fe . We 
obtain a = 0.091 ±0.015 and b = 0.68 ±0.07. These results are plotted in Fig. 
2. This provides a simple empirical expression relating the coupling constant 
H and the Tsallis entropic parameter q. The limit q — > corresponds to a 
strongly correlated liquid which is able to achieve a high degree of order at 
the vicinity of the macroion surface. The weak coupling limit, reproduced tak- 
ing q — > 1 corresponds to a disordered liquid lacking of this ability to achieve 
order. It is worthy noting that the effect of the increasing correlations leads 
to structures in which certain microstates are privileged. These structures can 
be rigorously described by BG statistics, although this is a challenging math- 

ul computational resources. 
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12j work accurately 



ematical problem in most cases requiring power 
Some of the methods existing in the literature 
in the strong or weak coupling regimes, which correspond, approximately, to 
each of the two plateaus in Fig. 2. In this paper we have provided an alter- 
native and unified method, which provides complementary insight, providing 
accurate expressions for the counterion concentration profiles valid for arbi- 
trary coupling. The tricky intermediate region is also accurately accounted 
for in the theory presented. Coupling is economically introduced through the 
parameter q which is related to the relevant constant controlling the strength 
of the pair interaction term in the Hamiltonian lol ]. 

It is worthy insisting on the fact that the exact many particle BG canoni- 
cal partition function is the appropriate tool for the problem addressed here. 
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Previous works and methods cited above are quite respectable and do con- 
stitute rigorous approaches that have provided many relevant insights in the 
physics involved, although restricted to certain coupling regimes. We cast no 
doubt on the validity and suitability of these methods. The issue addressed 
here is a complimentary view that provides additional insight, approximat- 
ing the exact many-particle BG distributions mean field single-particle Tsallis 
distributions. Concretely, closed analytical expressions for the counterion con- 
centration profiles valid for the whole spatial range and arbitrary coupling 
have been provided in excellent agreement with previous numerical results. 
Previously calculated analytical limits are also exactly reproduced by the the- 
ory. 

In the preceding treatment we have considered pointlike counterions. The 
theory can be applied to other geometries through the correct choice of the 
Laplacian in Eq. (12) and to a wide range of problems as, for example, poly- 
electrolyte multilayers. These questions, as well as studying charge inversion 
in the present framework, will be the aim of future work. 

We want to acknowledge an anonymous referee for helpful comments. Financial 
support from the MCYT (Ministry of Science and Technology of Spain) and 
FEDER under Project No.MAT2002-00646 is also gratefully acknowledged. 
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Figure Captions 



(1) Counterion concentration profiles in dimensionless variables. Symbols cor- 
respond to the results of Monte Carlo simulations from Ref . 10J . The lines 
have been calculated from Eq.(16) for the values of q indicated in the 
figure. The dotted line corresponds to the solution of the PB equation 
(q = 1). The inset shows, in logarithmic scale, curves for increasing value 
of the dimensionless position where other cases of moderate coupling have 
been added. 

(2) Comparison between data taken from Table 1 and a function q = 1/(1 + 
aE) b for a = 0.091 and b = 0.68. 
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Fig. 2. 
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E q 


0.1 


0.99 


1 


0.92 


10 


0.65 


100 


0.2 


10 4 


0.02 


10 5 


0.005 



Table 1 

Values of q corresponding to those of the coupling constant S considered in Ref. 
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